##########
###MAPS###
##########
library(sp)
library(RColorBrewer)
library(ggplot2)
library(maptools)
library(scales)
library(ggplot2)
library(doBy)
library(plyr)
library(foreign)

#set working directory
setwd("C:/Projects/Replication files/Data/")
main.dat <- read.csv("geol_dat_analysis_upd_pop3.csv")
###Collapse data by province
prov.dat <- summaryBy(event + conf + viol ~ prov, keep.names = T, FUN=c("sum"), data=main.dat)
prov.dat2 <- summaryBy(irrig_max + hyd_pop1950 ~ prov, keep.names = T, FUN=c("mean"), data=main.dat)
prov.dat.fin <- join(prov.dat, prov.dat2)
prov.dat.fin$id <- prov.dat.fin$prov
prov.dat.fin$prov <- NULL

####Import the shapefiles
#set new working directory
setwd("C:/Projects/Replication files/Data/KEN_adm_shp/")
ohsCol2 <- readShapeSpatial("KEN_adm1.shp")
ohsColI2 <- fortify(ohsCol2)
ohsColI2 <- data.frame(ohsColI2)

#set to old working directory
setwd("C:/Projects/Replication files/Data/")

##########################
###Total Conflict Event###
##########################
###Merge the by-province data
c1 <- prov.dat.fin[,c("id", "event")]
#Lag the province ids by one to fit the shapefile id framework
c1$id <- ifelse(c1$id==33, 32, 
                ifelse(c1$id==22,21,9))
#Join the two
ohsColI2.s1 <- join(ohsColI2, c1)
#ohsColI2.s1$event[is.na(ohsColI2.s1$event)] <- 0
###Plot comparative map
map_event <- ggplot() +
  geom_polygon(data=ohsColI2.s1, aes(x=long, y=lat, group = group,
                                     fill = event), colour ="black", size = 0.1) +
  scale_fill_gradient(low = "#FF9999", high = "#990000", space = "Lab", na.value = "white", guide = "colourbar",limits=c(0, 140)) +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1))+ # make title bold and add space
  labs(title = "All Political Violence Events By District", fill = "N. Events") +
  labs(x="",y="",title="") 
jpeg("map_event.jpeg", width = 5.5, height = 5, units = 'in', res = 500)
map_event
dev.off()

####################
###Irrigated Area###
####################
###Merge the by-province data
c3 <- prov.dat.fin[,c("id", "irrig_max")]
#Lag the province ids by one to fit the shapefile id framework
c3$id <- ifelse(c3$id==33, 32, 
                ifelse(c3$id==22,21,9))
#Join the two
ohsColI2.s3 <- join(ohsColI2, c3)
###Plot comparative map
map_irrig <- ggplot() +
  geom_polygon(data=ohsColI2.s3, aes(x=long, y=lat, group = group,
                                     fill = irrig_max), colour ="black", size = 0.1) +
  scale_fill_gradient(low = "#56B1F7", high = "#132B43", space = "Lab", na.value = "white", guide = "colourbar",limits=c(0, 670)) +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1))+ # make title bold and add space
  labs(title = "Average Irrigated Area (1950) By District", fill = "Sqrd. Km.") +
  labs(x="",y="",title="") 
jpeg("map_irrig.jpeg", width = 5.5, height = 5, units = 'in', res = 500)
map_irrig
dev.off()

################
###Population###
################
###Merge the by-province data
c4 <- prov.dat.fin[,c("id", "hyd_pop1950")]
c4$hyd_pop1950 <- c4$hyd_pop1950/1000
#Lag the province ids by one to fit the shapefile id framework
c4$id <- ifelse(c4$id==33, 32, 
                ifelse(c4$id==22,21,9))
#Join the two
ohsColI2.s4 <- join(ohsColI2, c4)
###Plot comparative map
map_pop <- ggplot() +
  geom_polygon(data=ohsColI2.s4, aes(x=long, y=lat, group = group,
                                     fill = hyd_pop1950), colour ="black", size = 0.1) +
  scale_fill_gradient(low = "#56B1F7", high = "#132B43", space = "Lab", na.value = "white", guide = "colourbar",limits=c(258, 924)) +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        plot.title = element_text(lineheight=.8, face="bold", vjust=1))+ # make title bold and add space
  labs(title = "Average Population, 1950 (Estimated, in Thousands)", fill = "People (1000s)") +
  labs(x="",y="",title="") 
jpeg("map_population50.jpeg", width = 5.5, height = 5, units = 'in', res = 500)
map_pop
dev.off()
